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Abstract 

The solutional origin of limestone caves was recognized over a century ago, 
but the short penetration length of an undersaturated solution made it seem 
impossible for long conduits to develop. This is contradicted by field ob- 
servations, where extended conduits, sometimes several kilometers long, are 
found in karst environments. However, a sharp drop in the dissolution rate 
of CaCOs near saturation provides a mechanism for much deeper penetra- 
tion of reactant. The notion of a "kinetic trigger" - a sudden change in rate 
constant over a narrow concentration range - has become a widely accepted 
paradigm in speleogenesis modeling. However, it is based on one-dimensional 
models for the fluid and solute transport inside the fracture, assuming that 
the dissolution front is planar in the direction perpendicular to the flow. 
Here we show that this assumption is incorrect; a planar dissolution front 
in an entirely uniform fracture is unstable to infinitesimal perturbations and 
inevitably breaks up into highly localized regions of dissolution. This pro- 
vides an alternative mechanism for cave formation, even in the absence of 
a kinetic trigger. Our results suggest that there is an inherent wavelength 
to the erosion pattern in dissolving fractures, which depends on the reaction 
rate and flow rate, but is independent of the initial roughness. In contrast 
to one-dimensional models, two-dimensional simulations indicate that there 
is only a weak dependence of the breakthrough time on kinetic order; local- 
ization of the flow tends to keep the undersaturation in the dissolution front 
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above the threshold for non-linear kinetics. 
Keywords: dissolution, speleogenesis, hydrology 



1. Introduction 



The notion that caves result from dissolution by an aqueous solution of 
carbon dioxide was already present in the works of Lyell and Thirria in the 
1830s (jShawl . 120001 ). but it was not until 100 years later that these ideas 
were developed mathematically A quantit ative model of the dissolution of a 
single limestone fracture was developed by IWeyll ( 119581 ). taking into account 
chemical kinetics and solute transport. His theory led to the paradoxical con- 
clusion that water flowing through a limestone fracture becomes saturated 
with calcium ions over length scales of the order of centimeters, so that lime- 
stone caves should not exist at all (IWhite and Longyearl . 119621 ) . A possible 
resolution of this paradox was proposed by I White! ( 119771 ). who noted that the 
existence of large cave systems may be explained by the sharp drop in the 
dissolution rate of CaC0 3 near saturation; this is frequently referred to as 
the "kinetic trigger" mechanism in the speleogenesis literature. Dissolution 
of blocks of calcite under laboratory conditions confirms that there is a rapid 
drop in reactio n rate near saturation, apparently because o f impuri t ies in 



natural calcite (Plum mer and Wiglevl . Il976t iDreybrodtl . Il990t iPalmerl . 11991 
Svensson and Drevbrodtl. Il9 92). 



Se yeral calcul a tions (IDreybrodtl . Il990l ; IPalmerl . Il99ll ; I Groves and Howard! 



19941 ; IDreybrodtl . Il996l ) have shown that the kinetic trigger hypothesis al- 
lows for the development of deep conduits, but they are all based on a one- 
dimensional model of fracture dissolution. The initial fissure is approximated 
by two parallel planes, and all the relevant variables (aperture, fluid volume 
flux, and solute concentration), depend only on the distance from the inlet 
(see Fig. |TJ). However real fractures are never one-dimensional; in this paper 
we show that even a tiny variability in fracture aperture inevitably leads to 
an instability at the reaction front, which then breaks up into highly local- 
ized regions of dissolution. Thus, the premise that a smooth fracture will 
open uniformly across its width, which underlies current models of speleoge- 
nesis, is faulty. A correct description of fracture dissolution must include the 
variation in aperture across the lateral direction as well. 
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h = h(x) 



Figure 1: Dissolution of a one-dimensional fracture; fluid flow is in the x direction and 
the fracture surfaces dissolve in the normal (z) direction. The fracture surfaces can be 
assumed to be located at ±/i/2. 



2. Breakthrough times in one and two dimensions 

In this paper we reevaluate the standard mathematical model for the early 
stages of cave formation, describing the dissolution of a calcite fracture by 
surface water draining through it to a lower-lying water table. We analyze the 
coupled equations for fluid flow, reactant transport and surface dissolution 
to show that the evolution of the fracture aperture is an inherently two- 
dimensional process, even when the initial aperture field is spatially uniform. 

2.1. Flow and transport: ID model 

In studies of fracture dissolution, and particularly in theoretical inves- 



frequently used (e.g. 


Drevbrodt. 


1996; Dijk and Berkowitz 


. 1998) 



1990; I Groves and Howard! . Il994t iDreybrodt 



Fig. [U the fracture aperture h(x, t) is assumed to depend on a single spatial 
variable, the distance from the inlet. The flow rate, q(t), is then independent 
of position, 



Ap 

W) 



r(t) = 12// 



dx 



h(x, ty 



(i) 



where Ap is the difference between the inlet and outlet pressures, L is the 
fracture length along the flow (x) direction, and r is the integrated flow 
resistivity. 
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The concentration of Ca 2+ ions in the fracture, c(x,t), is described by a 
convection-diffusion equation 



q^--^( hD xx ^-) = 2R(c), 
ax ax \ ax J 



(2) 



where D xx is the dispersion coefficient and R{c) is the reactive flux from 
the dissolving calcite. The factor of two in the dissolution rate comes from 
combining erosion at the upper and lower fracture surfaces. It is assumed that 
the inlet stream is pure water, c(0, t) = 0, and the outlet stream is a saturated 
solution, c(L,t) = c sat . Finally, the aperture evolution is determined from 
the local reactive flux 

cj± = 2R(c) (3) 
where c so \ is the molar concentration of the solid phase. 
2.2. Dissolution and breakthrough: ID model 



E arly theories of karstification assumed linear dissolution kinetics (IWeyll . 
19581 : IWhite and Longvearl . Il962h . 



R(c) = k(c sat - c) 



(4) 



In this case the undersaturation decays exponentially into the fracture, c sat — 



c ~ e 



-x/l f 



, with a penetration length, 



l p = q/2k. 



(5) 



Dissolution is restricted to a narrow region near the inlet, x ~ l p , which limits 
the growth of conduits. In fractured carbonate formation s l p is typically les s 
than a meter ( jWhite and Longvearl . Il962t lAtkinsonl . Il968t iDreybrodtl . Il990l ) . 
White! ( 1l977l ) proposed that field observations of extended conduits in karst 
environments, sometimes several kilometers long, might best be explained by 
a sharp drop in mineral dissolution rate as saturation is approached. This 
"kinetic trigger" mechanism can be modeled as a switch from linear to higher 
order kinetics at a threshold concentration c t h'- 



R{c) 



k{c S at c )? 
knifisat ~~ c )' 



C < C th 
C > C th - 



(6) 



A non-linear kinetic law (n > 1 in Eq. (jSJ)) results in an algebraic decay of 
the concentration profile at large distances from the inlet, c sat — c ~ a; 1 / 1- ™, 
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Figure 2: One-dimensional dissolution of a 20m fracture with a uniform initial aperture 
h(x) = ho = 0.2mm; Pe = 100 and Da = 5 x 10~ 4 . The volumetric flux q(t), scaled by 
the initial flux in the fracture qo, is shown for linear and non-linear kinetics. The inset 
shows the concentration profile from Eqs. @ and ©, with cth = 0.9c sa t. 



and the fracture opens over its entire length. Concentration profiles for linear 
and non-linear kinetics are illustrated in the inset to Fig. |2j In these (and 
subsequent) calculations we take typical parameter values for c alcite disso- 
lution : n = 4, k = 2.5 x 10~ 5 cms _1 and c sat = 2 x 10~ 6 Mcm~ 3 f Drevbrodtl . 



19961 ). The reaction rate fc 4 = k(c sat — c^) -3 is adjusted so that R(c) is con- 
tinuous at c = c t h, and the threshold concentration itself was set to 0.9c sa t. 
The molar concentration in the solid phase was taken as 0.027Mcm -3 , based 
on the mass density of pure calcite. 

In this paper we take the initial fracture aperture ho = 0.2mm, which 
means that the product of Peclet number, Pe = qo/D, and Damkohler num- 
ber, Da = kh /q , is about 0.05. Here q is the initial value of the volumetric 
flux in the fracture and the solute diffusion constant D ~ 10 _5 cm 2 s _1 . 

Figure [2] shows typical breakthrough curves for linear and non- linear ki- 
netics, obtained by solving Eqs. [THS] numerically (see Appendix A for details). 
The volumetric flux in the initial fracture is about 10 _3 cm 2 s _1 , corresponding 
to an hydraulic gradient of 0.01. The penetration length, L w 20cm, is much 
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less than the length of the fracture (L = 20m). Thus with linear kinetics the 
opening of the fracture aperture is slow, with the breakthrough time scaling 
exponentially with the length of the fracture. Breakthrough occurs much 
earlier in the non-linear case, because t he slower kinetics allow a deepe r pen- 



etration of reactant into the fracture (jWhitd . 119771 ; iDreybrodtl . I199CH ) . The 



time scale for dissolution can be characterized by the growth of the fracture 
aperture at the inlet, h(0,t) = h + 2ktc sat /c so i, which is independent of the 
kinetic model and the dimensionality of the fracture. We define 



r = h c so i/2kc sat 



(7) 



as the time it takes for the initial fracture aperture at the inlet to double. 
For an 0.2mm fracture, r = 5.4 x 10 6 s or approximately 0.17 years. 

2.3. Flow and transport: 2D model 

Although one-dimen sional models are simple and mathematically tractable, 



labor atory experiments ([Durham et al.l . l2001t iGouze et all l2003t iDetwiler et al 



20031 ) have shown that in most cases fracture dissolution is non-uniform in 
the direction transve rse to the flow, with highly localized two-dimensional 
dissolution patterns ( lHanna and Rajaraml . Il998[ ). The non-linear dynam- 
ics underlying th is more complicated behavior can be probed using depth - 
averaged models (jCheung and Rajaraml . 120021 ; IDetwiler and Rajaraml . 120071 ). 
in which the fluid velocity and reactant concentration are still averaged over 
the fracture aperture (z direction in Fig. [I]), but can vary in the lateral 
(y) direction. Coupled equations for the fluid volume flux q(x,y,t), depth- 
averaged concentration of dissolved solids c(x,y,t), and aperture h(x,y,t) 
are solved simultaneously: 



Vq = 0, q 



h 3 Vp 

qV ■ c - V • (hD ■ Vc) = 2R(c) 
CsoAh = 2R(c), 



(8) 



where D(h) is the solute dispersion tensor. Details of our numerical solution 
of Eq. can be found in Appendix A 



In Sec. [3] we will show mathematically that the one-dimensional dissolu- 
tion profiles described by Eqs. flTJ)-(J3]), are unstable. Uniform concentration 
profiles inevitably break down into more complex dissolution patterns, which 
can be described by Eq. (jSJ). Simulations based on Eq. (jBJ) have proven to be 
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Figure 3: Dissolution of a smooth fracture 20m x 10m x 0.2mm; the parameters Pe = 100 
and Da = 5 x 10~ 4 . One-dimensional simulations of the volumetric flux q(t) are compared 
with two-dimensional simulations of the average flux q a vg{t) — W^ 1 f q x {x, y, t)dt. The 
fluxes are scaled by the initial flux in the fracture q$. 



succes sful in predicting the large sc ale dissolution patterns in artificial frac- 
tures (IDetwiler and Raj araml . 120071 ). although at small scales, comparable to 
the correlation length in the aperture field, the dissolution is more uniform 
than the experimental observations. For such fine details three-dimensional 
simulations are necessary, including the variation in flow and concentration 
fields across the aperture (jSzymczak and Laddl . 120091 ) . 



2-4- Dissolution and breakthrough: 2D model 

We repeated the one-dimensional simulations shown in Fig. [2] using a 
two-dimensional fracture, 20m x 10m, superimposing a small random aper- 
ture field, with a variance a = 20nm, over the original uniform aperture 
(0.2mm). Time-dependent flow rates from ID and 2D simulations are com- 
pared in Fig. [31 Two-dimensional simulations show a significant reduction 
in the breakthr ough time in comparison w i th ID models, in agreement with 
earlier results (IHanna and Raj araml . Il998l ; ISzymczak and Laddl . 120091 ) . In- 
terestingly, the added spatial dimension has a larger effect on the break- 
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Figure 4: Concentration profiles in a dissolving fracture. Contour plots of the normalized 
undersaturation, (c sat — c)/c sa t, are plotted at successive times (from left to right); t = 
25t, t = 50t, t = 75t, and t — lOOr. The flow direction is from left to right; the 
Peclet and Damkohler numbers are Pe = 100 and Da = 5 x 10~ 4 . A video of the 
complete time sequence is included as an ancillary file to this article and also accessible 
at http : //www . f uw . edu . pl/~piotrek/ c . mpg. 



through time than the kinetic trigger. In the one-dimensional simulations, 
non- linear kinetics reduces the breakthrough time by an order of magnitude, 
from 2800 years to 270 years. However, in the two-dimensional simulations, 
the breakthrough time is less than 30 years, regardless of the kinetics. In 
fact the non-linear kinetics makes only a small difference here, reducing the 
breakthrough time from about 28 years (linear) to 26 years (non linear). 

The reason why dimensionality is crucial can be seen in the concentra- 
tion maps shown in Fig. HJ The initially planar dissolution front breaks down 
into increasingly focused regions of dissolution, which rapidly advance into 
the fracture, causing breakthr ough at much earlier times than in the homoge- 
neou s (one-dimensional) case (jHanna and Rajaraml . ll998l ; ISzymczak and Ladd . 



20091 ) . Flow focusing keeps the concentration in the advancing front high, 



above the threshold for non-linear kinetics; this explains the relatively small 
difference in breakthrough times in the two-dimensional simulations. 

In the course of extensive numerical investigations of fracture dissolution, 
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we noticed that a planar dissolution front was always unstable if the fracture 
length, L, was sufficiently large in comparison with the penetration length, 
l p . Our simulations suggested that on geophysical scales, where L ^> l p 
dissolution will be inherently unstable. Mathematical analysis, described 
in|Appendix C] confirms that fracture dissolution is always two-dimensional, 



and cannot be usefully approximated by one-dimensional models. 

3. Instability of dissolution in a uniform aperture fracture 

The stability of a uniform dissolution front can be investigated from a 
minimal model, assuming linear kinetics, and neglecting the dispersive term 
in Eq. ©: 

V -q = 0, Q=-jk — , 

12 t* (a) 
q-Vc = 2k(c sat - c) {y) 

c so id t h = 2k(c sat - c). 

In Sec. 12.41 we showed (Fig. [3]) that once a two-dimensional model of frac- 
ture dissolution is adopted, variations in kinetic order are of only quantita- 
tive rather than qualitative importance to the predicted breakthrough time. 
Moreover, on the scale of penetration length, the ratio of the diffusive flux, 
Dhc/lp, to the convective flux, qc/l p , is of the order of Pe _1 Da. This param- 
eter is small in typical limestone fractures, as shown by a consideration of 



fracture apertures and hydraulic gradients ( |Appendix B[ ); thus Eq. ([9]) con- 
tains all the essential physics of the instability. However, in gypsum karst 
the reaction rate is much higher ( Da > 1) and the penetration length equiva- 
lently smaller; in this case diffusive flux is significant and should be included, 
as in Eq. (ICTTil) . 

In the initial stages of dissolution, when the penetration length is much 
less than the length of the fracture, l p <C L, the one-dimensional flow and 
concentration fields are essentially time independent, 

9(*)«ft = -^. 4x,t)ac {x) = Crt{l-e-' l >). (10) 

The pressure gradient p' = d x p is constant far from the inlet and q is inde- 
pendent of x by continuity; the subscript is used to denote initial values. 
The aperture grows linearly in time, with a rate that decays exponentially 
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with the distance from the inlet, 



h{x,t) = h [l + ). (11) 



The flow rate and concentration field remain constant in time until the pene- 
tration is deep enough to perturb the far field pressure gradient and increase 
the flow rate in the fracture. 

However, an initially smooth fracture does not evolve according to Eq. (fTTj) : 
in fact the one- dimensional solution is exponentially unstable to perturba- 
tions in h, 5h(x, y, t) ~ e ut , which rapidly outrun the uniform front suggested 
by Eq. ( ITT]) . A linear stability analysis of Eqs. §§§ leads to a dispersion re- 
lation for the growth rate u as a function of the wavelength of the pertur- 
bation A, as shown in Fig. |5j The particulars of the analysis can be found 
in |Appendix C| Here we just mention one subtle but important detail; the 



base state, Eq. (ITT]), is itself time-dependent. Th e stability of nonauto n omou s 



systems is in general a very difficult problem ( iFarrell and Ioannoul. Il996l); 
here however we follow a relatively simple, approximate approach (jTan and Homsy 



19861 ) in which the base state is frozen at a particular time and the growth rate 
is determined as if the base state were time-independent (the quasisteady- 
state approximation). The dispersion curve in Fig. [5] was obtained by freezing 
the base state at to = 0. We will examine this approximation in more detail 
in Sec. SJ 

The most striking result of the linear stability analysis is that there is 
a maximal growth rate, u max = 0.79r _1 , at a wavelength X max = 4.74/ p 
(Fig. [5]). An individual fracture can therefore develop a strongly heteroge- 
neous permeability during dissolution, with an inherent length scale that 
depends on the kinetics and flow rate (via l p ), but not the initial topogra- 
phy. Furthermore, there is no lower limit to the reaction rate for unstable 
dissolution. Only the wavelength and penetration depth are affected, scaling 
with the inverse of the Damkohler number. 

In laboratory experiments it is frequently the case that the sample length 



is less than l p , in which case uniform dissolution is observed ( jFredd and Foglerl . 



19981 ) . But on scales of geophysical relevance, dissolution of carbonate rocks 
will always be unstable; calcite for example would have characteristic wave- 
lengths Xmax ~ lm, while in dolomite X max ~ 10m. It is worth noting that 
the instability in the dissolution front continues to grow in systems confined 
to widths smaller than l p , although at a reduced rate. Figure |5] shows that 
there is a long tail to the dispersion relation at short wavelengths, with a 



10 




Figure 5: Growth rate of the instability in a smooth fracture. The time scale r = h n /2kj. 
The inset figure shows the growth rate over a larger range of wavelengths, including the 
effects of lateral diffusion as measured by the product Pe^ 1 Da. 
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weak power-law dependence uj ~ A. Thus in narrow fractures, where the 
width W < X max , a single channel develops which contains all the flow in the 
fracture. 



The instability analysis in Appendix C can be generalized to include lat 



eral diffusion (in the y direction), which is expected to play a more important 
role in the dynamics than axial diffusion because axial transport is usually 
convection dominated. On scales comparable to the penetration length, the 
relative magnitude of diffusive and convective fluxes is given by the param- 
eter Pe _1 Da, as can be seen in Eq. (1C.14|) . The inset to Fig. [5] shows the 



dispersion relation for several different values of the product Pe _1 Da. Lat- 
eral diffusion reduces the growth rate of short wavelength perturbations and 
shifts Xmax towards longer wavelengths. However, it does not prevent the in- 
stability developing and the growth rate remains positive for all wavelengths. 

4. Numerical simulations of the instability 

We have confirmed the key predictions of the instability analysis by nu- 
merically solving Eqs. ([H]) with linear reaction kinetics. In comparison with 
the minimal model considered in the previous section, we have here in- 
cluded solute dispersion and the effect of diffusion on the dissolution kinetics, 
Eq. (1A.5|) . Beginning from a uniform fracture (ho = 0.2mm), with a small 



random roughness superimposed (a = 10~ 4 h ), we see a single sinusoidal 
mode developing in the dissolution front, as illustrated in Fig. |H The wave- 
length and front thickness are inversely proportional to Damkohler number, 
as shown in Fig. [61 confirming that penetration length, l p = /io/2Da, is the 
only important length scale in the early stages of dissolution. Concentration 
profiles, such as those illustrated in Fig. |6l can be Fourier transformed in the 
lateral (y) direction to determine the wavelength of the instability as a func- 
tion of the Damkohler number. An analysis of the number of wavelengths, 
W/X max , in a system of width W = 10m confirms that the fastest growing 
modes have a wavelength close to the theoretical prediction X max = 4.74/ p 
(l P — ho/2 Da), as can be seen in Fig. [3 

The growth rates of the instability can be obtained by analyzing fluctu- 
ations in concentration field, 



A c (x, t) = /5>(x, y, t) - c(x, t)] 2 . (12) 
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Figure 6: Concentration profiles of a dissolving fracture. Contour plots of the normalized 
undersaturation, (c sa t — c)/c so t, are plotted for different Damkohlcr numbers (from left 
to right); Da = 1.25 x 1(T 4 , 2.5 x 1CT 4 , 5.0 x 10~ 4 , and Da = 10~ 3 . The product of 
Peclet and Damkohler numbers is 0.05 in each case. The flow is from left to right and the 
time t = 50r. The variance in the aperture was reduced to 10~ 5 ho so that the sinusoidal 
perturbations in the front are more easily visible at larger Da. 
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Figure 7: Wavelength dependence of the dominant mode as a function of Damkohler 
number, Da. The solid circles show the number of wavelengths in a system of fixed width 
W = 10m, determined by a Fourier analysis of the concentration profile. The solid line is 
the theoretical prediction based on wavelength of the fastest-growing mode, X m ax — 4.74i p 
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Figure 8: Time-dependent fluctuations in concentration for different values of the 
Damkohler number, Da, and roughness, cr/ho; a is the variance in the aperture. The 
solid line shows the growth of the mode with an initial wavelength of A max , based on a 
numerical integration of the linearized equation, Eq. (|C.7j) 



We take the maximum value of A c (x,t) along the fracture, A c (i), as a mea- 
sure of the amplitude of the perturbation at time t. Results for A c (t), 
presented in Fig. [8j confirm that the instability is initiated as soon as the 
undersaturated solution enters the fracture, and that the growth is indeed 
approximately exponential in time. In reduced time units (t/r) the growth 
rate is independent of Damkohler number and roughness; only the amplitude 
of the fluctuations is affected by the initial variability in the fracture aper- 
ture. For comparison, we also plot the solution of the linearized initial value 
problem, Eq. ( 1C.7I) . which includes the time dependence of the base aperture 
field, Eq. ( TTTT) . We take the eigenf unction, Eq. ( 1C.12I) . corresponding to the 
fastest-growing eigenmode, u = 27il p /\ max , as the initial condition. Initially, 
fluctuations in a random aperture field grow more slowly than the fastest- 
growing mode, but after a while (~ lOr) a single mode is dominant, and the 
growth rates become similar. 

The linear stability analysis summarized in Fig. predicts an exponential 
growth of the instability. On the other hand, the numerical results in Fig. [H] 
indicate that the growth is less than exponential. This is because the time 
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Figure 9: Growth rate of the instability in a smooth fracture with the base state frozen at 
different times: to — (solid line), to — r (dashed line), and to = lOr (dotted line). The 
inset figure shows the growth rate of the mode with an initial wavelength of \ ma x — 4.74Zp, 
for different times: the linearized solution (solid line), including the time-dependence of 
the base state, Eq (IC.7p . is compared with the eigensolution analogous to Eq. (IC.13|) 
(dashed line) 

dependence of the base state weakens the instability as the fracture inlet 
opens up. Figure [9] shows the dispersion relation obtained from the quasi- 
steady state approximation, Eq. (IC.9j) . for different base states, frozen at 
t = (as in Fig. [5]), t = r, and t = 10r; the growth rate is reduced 
across the whole spectrum of wavelengths as time goes on. However, the 
peak growth rate remains at almost the same wavelength, independent of 
time, so a single mode (A ~ \ ma x) predominates until the onset of non-linear 
growth. 

The quasi-steady-state approximation, Eq. (1C.9I) . used to calculate the 
dispersion curves in Fig. M (and Fig. \5§ overestimates the growth rate when 
compared to the exact linearization in Eq. (10.7[) . The solid line in the inset 
to Fig. [9] shows the growth rate of the mode with an initial wavelength of 
^max as a function of time; in essence the time derivative of the solid line in 
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Fig. EJ The dashed line is the growth rate for the same wavelength (X max ) 
obtained by freezing the base state at successive moments of time. The 
discrepancy is about 40 % at to = and decrea ses thereafter, which agrees 



with the observation of iTan and Homsyl (119861 ) that the quasi-steady-state 
approximation becomes exact in the long-time limit. 



5. Discussion 

The dissolution of an initially smooth fracture is one dimensional at early 
times, but soon tiny random variations in fracture aperture lead to the devel- 
opment of small perturbations in the front. The strong wavelength selection 
(Fig. [5]) means that after a fairly short time, approximately 1 — 2 years in our 
model calcite fracture, a single mode becomes dominant, with a wavelength 
that is about 5 times the penetration length (Fig. [7]). As the perturba- 
tion becomes stronger, distinct channels emerge from the front as shown in 
Fig. HI The Ivideol illustrates this process dynamically; the shorte r channels 
are drained of reactant and cease to grow (IFredd and Foglerl . Il998l ). Thus the 
spacing between long conduits is not determined solely by the initial wave- 
length of the instability but also by the non-linear dynamics of interacting 
channels. Eventually all the flow is localized within a few active channels. In 
this paper we have shown that this is the expected behavior on sufficiently 
large scales. 

The notion of a kinetic trigger is an attractive explanation for the exis- 
tence of long conduits in calcite formations, where there is an abundance of 
evidence for a t ransition from linear to non- line a r kine tics (jPlummer and Wigleyi . 
19761 ; iPalmerl . Il991t ISvensson and Dreybrodtl . Il992l ). However, in gypsum 
there is good reason to think that the dissolution rate is roughly l inear 
in the undersaturation (jColombani and Bertl . 120071 ; IColombanil . 120081 ). al- 
though it has been suggested t hat there is a chang eover to non- linear kinet- 



ics in natural gypsum as well ( IJeschke et al.l . l200ll ). One-dimensional mod 



els of gypsum dissolution do not develop long con duits for typical fracture 
apertures and flow rates (IRaines and Dewerj . 1 19971 ), unless a kinetic trigger 



with a high-order (n = 4) rate law is again invoked (IDreybrodt et al.l . 12002 



Romanov et al. . 20031 ). 



The penetration length, Eq. (J5J), increases with time due to the increasing 
flow rate, which introduces a feedback mechanism that produces the charac- 
teristic upturn in q(t) near breakthrough (Figs [2] and [3]). In one-dimensional 
models an increase in flow rate can only arise from the growth of the mean 
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aperture. However in two- dimensions flow focusing provides a mechanism for 
a localized increase in flow rate; longer channels drain the flow from neigh- 
boring regions of the fracture, promoting their own growth at the ex pense 
of their neighbors (IFredd and Foglerl . ll998l ; lHanna and Rajaraml . Il998l ). The 
process can be un derstood by consider i ng th e pressure fields in channels of 
different lengths (jSzymczak and Laddl . 120061 ). Given enough time, only a 
single channel remains and at that point the flow focusing mechanism stalls; 
the flow rate in the remaining channel can now only grow by extending its 
length. Thus there is a maximum flow rate that a single channel can acquire 
by flow focusing, q max ~ q W/w, where w is the width of a single dissolution 
channel or wormhole. These observations suggest that the aspect ratio of 
the fracture and its orientation with respect to the flow may also play an im- 
portant role in determining the breakthrough time. If W > L the drainage 
basin is large and breakthrough will be rapid (as in the video), but if W <C L 
the flow focusing may stall with the dissolution front far from the outlet; in 
this case breakthrough will be dependent on reaction kinetics. 



D iscrete fracture networks are w idely used in reservoir modeling (ILong et al. 



1982 ; ISiemers and Drevbrodtl . Il998l ) and in the assessment of dam safety (IDreybrodt et al. 



20021 ; Romanov et al. . 20031 ) . because they provide an effective means to de- 



scribe localized regions of high permeability. However, each fracture within 
the network is modeled at the one-dimensional level, and the possibility of 
spatially localized dissolution within an individual fracture in the network 
is not taken into account. In fractured carbonates an instability in the dis- 
solution front can reduce the breakthrough time in an individual fracture 
by orders of magnitude, thereby completely changing the hydraulic prop- 
erties of the global fracture network. It is therefore important to develop 
a model of evolving fracture permeability that includes the inherent het- 
erogeneity of the dissolution process. Although the initial wavelength of 
the developing instability is controlled by the reaction kinetics and flow 
rate, interactions between the developing channels play the dominant role 
in the long-time evolution of the fracture aperture field. Water through- 
out most of the fracture is saturated with calcium ions; only in the narrow 
regions formed by the active channels is there any significant undersatu- 
ration. Thus, an adequate understanding of how single channels advance 
(Fig. [Hand http://www.fuw.edu.pl/~piotrek/cmpgvideo) would open up 
the possibility of new field-scale models of fracture dissolution, describing 
the growth of individual channels coupled together by long-range pressure 
fields. This could make it feasible to include heterogeneous dissolution of 
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individual fractures in reservoir modeling, since the pressure fields around 
long thin channels can be rapidly calc ulated using techniques of conformal 
mapping (IGubiec and Szymczakl . 120081 ) . 

Our analysis di f fers f r om previous work on instabilities i n por ous me- 



dia (jChadam et all Il986t ISherwoodl . Il987t iHinch and Bhattl . Il990h in that 



we take a time-dependent aperture as the base state, rather than a propagat- 
ing one-dimensional dissolution front. In a porous rock all the soluble ma- 
terial eventually dissolves and therefore a time-independent porosity profile 
can steadily advance into the material. However, in a fracture the aperture 
profile evolves along the whole length, with an effectively unlimited increase 
in aperture near the inlet. 



6. Conclusions 



The main conclusion from this work is that the dissolution of a single frac- 
ture is inherently two-dimensiona l; the one-dimens ional solutions frequently 
used in models of cave formation (jDreybrodtl . Il996l ) are unstable to infinites- 
imal perturbations. This means that the kinetic trigger mechanism is not 
a prerequisite for the development of long conduits. Although the develop- 
ment of instabilitie s in fracture dissolut i on has been demonstrated by nu 



merical simulation (IHanna and Rajaraml . Il998l ; I Cheung and Raj ar ami . 12002 



Szymczak and Laddl . 120061 . [20091) , this is the first time, to our knowledge, 



that the equations for fracture dissolution have been shown mathematically 
to have unstable solutions. We have further shown that the instability devel- 
ops with a well-defined wavelength, which depends on reaction kinetics and 
flow rate but is insensitive to the initial roughness of the fracture. Subse- 
quently, the localized regions of dissolution advance in the fracture in ways 
that may eventually be understood in terms of modified models of Laplacian 
growth. 
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Appendix A. Simulation method 



The depth-averaged equations, Eq. (jSj), were solved on a uniform grid 
with a resolution of 1cm; further grid refinement did not affect the results 
significantly. Periodic boundary conditions were imposed in the lateral (y) 
direction, while the inlet and outlet conditions were the same as in the one- 
dimensional case, Sec. 12.11 We included Taylor dispersion in the axial diffu- 
sion, D xx = D + ql/210D, but its effect on the erosion rate was negligible; 
in the lateral direction we took D yy = D. One-dimensional calculations used 
the same code, but with a single grid point in the lateral (y) direction. 

Calcite dissolution kinetics are more properly expressed in terms of the 
calcium ion concentration at the fracture surface, c s , rather than the bulk 
concentration as in Eq. (J6]): 

R( Cs ) = I k[ ° sat ~ Cs \ n Cs < Cth (A.l) 

The question then is how to relate the surface concentration, c s , to the bulk 
concentration, c; for this we must consider how diffusion limits the rate of 
mass transfer from the fracture surface. Diffusive transport of reactant across 
the fracture aperture can be ap proximated using an effective mass transfer 



coefficient or Sherwood number ( jBird et al.l . 1200 ll ). 



„ D Sh . , . . 

Rdiff = ^ fr (c s -c). (A.2) 

The Sherwoo d number, Sh, depends on reaction rate but the variation isrela- 



tively small (IHayes and Kolaczkowski I1994J ; iGupta and Balakotaiahl . 120011 ) , 



bounded by two asymptotic limits: constant flux, where Sh = 8.24, and 
constant concentration, where Sh = 7.54. In the numerical calculations we 
used the approximate value Sh = 8. 

In sufficiently narrow apertures dissolution kinetics are reaction limited, 
R <C Rdiff, and c s ~ c which gives the reactive flux in Eq. (@|. However, as 
the fracture opens the reaction rate becomes hindered by diffusive transport 
of reactant across the aperture. Finally, for kh/DSh ^> 1, the dissolution 
can become entirely diffusion limited, Rdiff <^ R- In that case the surface 
concentration approaches that of a saturated solution c s ~ c sat and 

.DSh . , . 

Rdiff = ~^( c ^t ~ c). (A.3) 
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Thus the reactive flux is again of the form of Eq. (J3J) but with an effective rate 
constant ko = D Sh/2/t. Weyll ( 1958 ) assumed diffusion-limited kinetics, but 



for a typical fracture aperture (10 _2 cm) kp is 1 — 2 orders of magnitude larger 
than k. In fact reaction-limited kinetics persist, with small modifications 
from the effects of diffusion, up until fracture apertures of the order of 1cm. 

An effective reaction rate, including the effects of diffusion, ca n be ob- 

tained by balancing the reactive flux, R, with the diffusive flux, Rdiff (jSzymczak and Ladd 



200% : 

^( Cs -c)=R(c s ). (A.4) 
In the case of linear kinetics we then obtain an effective rate constant 



k 

1 + {2kh/D Sh) ' 



k eff = i ; TTTTTTTToiIT' ( A - 5 ) 



which is valid for all apertures. For more reactive materials, such as gypsum, 
transport corrections are important, even for narrow apertures. In some cases 
the kinetics is entirely transport limited, as in Eq. (I A. 3 1) 

In the most general case, non-linear kinetics combined with diffusive mass 
transfer, we can combine Eqs. flA.lj) and (IA.4j) to obtain the reactive flux in 



terms of the bulk concentration c. Defining the dimensionless quantities: 

Csat c c sa t c s 2,kh . . 

C = , C s = , a = ——, (A.6) 

c sat ~ c th C sat ~ C t h U Oil 



Eq. flA.4j) becomes 

aC n s + C s - C = 0, (A.7) 

where C is known and C s is to be found. The numerical simulations use 
Eq. (IA.7p . but the instability analysis (Sec. E]) assumes reaction-limited ki- 
netics, Eq. 06]). The difference is small during the initial stages of calcite 
dissolution. 

We used a flux-conserving discre tization to solve for the pressure field, 



which avoids "numerical saturation" (jGroves and Howardl . ll994tlHanna and Rajaram , 



19981 ). even when the grid is relatively coarse. The scalar fields, aperture, 
pressure, and concentration are defined at the nodal points and gradients 
of the scalar fields are then naturally calculated on the edges of the cells 
surrounding each node. The divergence of the flux is then automatically cal- 
culated at the nodal position. However, the convective flux of reactant was 
calculated by upwind differencing to ensure stability. 



21 



The li near equations f or the pressu re field were solved using the MUMPS 



package (lAmestoy et all l200ll . 120061 ). and the fluid volume flux was then 



calculated by differencing. For linear reaction kinetics, the concentration field 
can be solved directly, but for non-linear kinetics Newton-Raphson iteration 
was employed, using MUMPS to solve the linear system at each step. Within 
each iteration of the bulk concentration field, the surface concentrations, 
Eq. ( lA.7p . must be determined iteratively as well. The transition between 
linear (C s > 1) and non-linear (C a < 1) kinetics corresponds to C = 1 + 
a, and thus the appropriate branch of Eq. (1A.1I) . can be identified from 
the bulk concentration. Once the concentration field has converged, the 
aperture field is updated based on the local erosion flux, 2R(c s ). We used 
an explicit midpoint method to determine the time evolution of the aperture 
field, requiring two complete cycles of the concentration solver at each time 
step. 



Appendix B. Transport and reaction parameters in typical frac- 
tures 



F r acture apertures are between . 005cm and 0.1cm (IMotyka and Wilkl . 
1984] ; IPaillet et all Il987l : IPrevbrodtl. Il996h. and hydraulic gradients are of 
the order of 10" 3 to IP" 1 (IPalmerl . Il99ll bilk and Berkowitzl . Il998t ). This 
gives a range of characteristic flow velocities in undissolved fractures from 
10 _4 cms _1 to lcms -1 . The corresponding Peclet numbers are 0.05 < Pe < 



10 4 , taking the solute diffusion coefficient D ■ 
rate f or limestone is u sually in the range 10" 



10 



9 - 
cm s 



5 cms 1 — 10 



19911 : iDrevbrodt 



l . The dis solution 
4 cms _1 ( IPalmerl . 



19961 ) . which leads to the Damkohler numbers in the range 
10 -5 < Da < 1. Thus in most calcite fractures the parameter Pe -1 Da is 
small and the diffusive flux can be neglected, as in Sec. [3j Diffusion plays a 
more prominent role in dissolution of fractured gypsum , where the reaction 
rates are of the order of O.Olcms -1 ( IJeschke et all 120011 ). and thus Pe -1 Da 
is 2 — 3 orders of magnitude larger than in limestone systems. 



Appendix C. Linear stability analysis 

Equation can be non-dimensionalized by scaling length by the pene- 
tration length, x — > 2kx/qo, and time by the time taken to double the initial 
inlet aperture, t — > 2ktc sat /hoc so i. It will be convenient to rewrite the con- 
centration field as a normalized undersaturation, c — > {c sa t — c)/c sa t- Then, 
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scaling the remaining variables (h, q, Vp) by their (constant) initial values 
(/t , Qo, p'o), we can rewrite Eq. ((H]) in dimensionless form, 

V • q = 0, <? = h 3 Wp, , Q ^ 

qr • Vc = — c, <9 t /i = c. 

The one-dimensional solution of (IC.lj) is taken as the base state: h b (x,t) = 
l+te~ x , q b = x, and c b = e~ x , where x indicates a unit vector pointing in the 
x direction. We consider perturbations about the base state, c = c b + 5c; h = 
hf, + 5h etc. and keep just the linear terms: 

V ■ 5q = 0, 5q = 3h b 1 5hx + hfVSp, , , 

fiqxd x Cb + d x 5c = —5c, d t 5h = 5c. 

The pressure variation can be eliminated from the equations for the vol- 
ume flux by cross differentiation, 

d y 5q x - d x 5q y = 3h^ (d y 5h - h' b 5q y ) . (C.3) 

where h' b = d x hb is the spatial derivative of the base aperture field. Com- 
bining this equation with the incompressibility equation, d y 5q y = —d x 5q x 
eliminates 5q y as well, 

(d 2 x + d 2 y )5q x = 3h b l {d 2 y 5h + h' b d x 5q x ) . (C.4) 

Finally, rewriting the transport equation from (10.2|) in the form 5q x = 
d x (e x 5c) and substituting into (10.4[) leads to equations for the linearized 
variations in the concentration and aperture fields, 

[h b (d 2 x + d 2 y ) - 3h' b d x ] d x (e x 5c) = W 2 5h, 3 t 5h = 5c. (C.5) 

We now examine the stability of a periodic perturbation in the aperture 
and concentration fields; 

5h = g(x, t)e~ x sm{uy), 5c = d t g(x, t)e~ x sm(uy), (C.6) 

where the (dimensionless) wavevector u = 2ttI p / A and A is the wavelength of 
the perturbation. This ansatz satisfies Eq. (IC.5I) provided that 

[h b (d 2 x - u 2 ) - 3h' b d x ] d t d x g + 3u 2 c~ x g = 0. (C.7) 
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In addition to the boundary condition g(Q, t) = 0, there are boundary condi- 
tions from the uniformity of the flow across the inlet and outlet, 5q y (0, y, t) = 
Sq x (x -> oo, y,t) = 0, or 



d 2 x g(0,t) = 0, d x g(x oo,*) = 0. 



(C.8) 



Equatio n (IC.7D does not have an eigensolution because of the time-dependent 
base state (jTan and Homsyl . Il986f ). but a numerical solution of Eq. flC.7j) is 
shown as the solid line in the inset to Fig. [9j Further analytic progress can be 
made by freezing the base state at a particular instant of time, to, replacing 
h b (x,t) by h b (t ) = 1 + t e~ x where £ is then kept constant: 



[h{toM - u 2 ) - 3h' b {t )d x ] d t d x g + 3u 2 e- x g = 0. 
In this case an eigensolution can be found of the form 

g(x,t) = g(x)e ujt , 
which gives an ordinary differential equation for g, 

[h b (t )(d 2 x - u 2 ) - 3ti b (t )d x ] d x g + fa-We-'g = 0. 



(C.9) 



(CIO) 



(C.ll) 



Here we present the solution for the simplest case, t = 0, for which h b (0) = 1 
and h' b (0) = 0. A similar but more complex solution can be obtained for 
arbitrary to an d results are shown in Fig. [9j 

The general solution of Eq. ( IC.llj) with to = is 



g(x) = A F 2 (1 + u, 1 - u; 3uo~ l u 2 e- x ) 

+ Be ux F 2 (1 + u, 1 - 2u; 3u 



l u 2 e~ x 



) 



+ Ce~ ux F 2 (l + u,l + 2u; 3u~ Ve") , (C. 12) 

where A, B, and C are constants and oF 2 (p, q; z) is a generalized hypergeo- 
metric function. The far field boundary condition d x g(x oo) = requires 
that B = 0, while the condition <?(0) = is then sufficient to determine the 
function g(x) to within an arbitrary constant, which is the initial amplitude 
of the perturbation. Imposing the final boundary condition, d%g(0) = 0, 



24 



gives a dispersion relation for u(u), 



u 2 F 2 (1 + u, 1 + 2m; 3uT V) + 

3(1 + 2u)u F 2 (2 + u, 2 + 2u; 3urV) 
+9w 2 F 2 (3 + u, 3 + 2m; 3w~V) ^2 (l + m, 1 - w; 3u;~V) = 
w 0^2 (2 + tt, 2 - w; So;" 1 ?/ 2 ) + 
3m 2 F 2 (3 + w, 3 - u; 3u~ V) ^2 (l + «, 1 + 2m; 3uT V) , (C.13) 

where oF 2 (p, q; z) = oF 2 (p, g; z) /T(p)T(q). The maximum growth rate (largest 
positive root) at each u from Eq. ( 1C.13j) is plotted in Fig. [5j 

The derivation of the dispersion relation can be generalized to include 
lateral diffusion (in the y direction). The convection-diffusion equation, in- 
cluding lateral diffusion and scaled in the same way as Eq. (10.2j) . contains a 
single parameter, the product Pe™ 1 Da = Dkho/q 2 ,, 

d x {e x 5c) - (2 Pe" 1 Da)d 2 y {e x 5c) = 5q x . (C.14) 

The analogous linear stability analysis leads to the generalization of Eq. flC.111) 
for finite diffusion. Simplifying to the case to = 0, 

(<9 2 - u 2 ) [8 X + (2 Pe" 1 Da)« 2 ] g + 3uj- l u 2 e~ x g = 0. (C.15) 

which can again be solved in terms of hypergeometric functions. The inset to 
Fig. [5] shows the dispersion relation for several different values of the product 
Pe" 1 Da. 

At long wavelengths the growth rate is linear in the wavevector, u — > 3u, 
while for small wavelengths an asymptotic analysis, including lateral diffu- 
sion, gives to — > 3PeDa _1 w~ 2 . However, the convective limit (no diffusion) 
appears to be singular, with oj ~ u~ x at large u. 
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